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Abstract. This paper presents a new computer code to solve the general relativistic 
magnetohydrodynamics (GRMHD) equations using distributed parallel adaptive mesh 
refinement (AMR). The fluid equations are solved using a finite difference Convex ENO 
method (CENO) in 3 + 1 dimensions, and the AMR is Berger-Oliger. Hyperbolic divergence 
cleaning is used to control the V ■ B = constraint. We present results from three flat space 
tests, and examine the accretion of a fluid onto a Schwarzschild black hole, reproducing the 
Michel solution. The AMR simulations substantially improve performance while reproducing 
the resolution equivalent unigrid simulation results. Finally, we discuss strong scaling results 
for parallel unigrid and AMR runs. 



1. Introduction 

The interaction of gravitational and electromagnetic fields together with rotation is believed 
to power the central engines of many astrophysical phenomena including relativistic jets in 
active galactic nuclei (AGN), other forms of black hole accretion, gamma-ray bursts (GRB), 
and core collapse supernovae. In addition, interactions between strong gravitational and 
electromagnetic fields are believed to result in the transport of angular momentum in accretion 
disks via the magnetorotational instability (MRI) and in the extraction of black hole energy 
via the Blandford-Znajek mechanism. 

The expected ubiquity of magnetic fields in the vicinity of strongly gravitating compact 
objects has spurred increased theoretical efforts to understand these astrophysical systems. 
The difficulties, however, are formidable as the physical laws describing these phenomena are 
nonlinear, evolutionary, and, in general, without simplifying symmetries. As a consequence, 
numerical simulation of these systems becomes crucial for better understanding them. Even 
then, there are likely considerable aspects of the physics that remain out of reach of current 
computational resources. For instance, it is not unrealistic to imagine that dissipative and 
radiative aspects of these problems will be important for accurately modeling certain types 
of phenomena. However, with the hope of capturing some of the relevant physics, a number 
of groups have begun developing methods and codes for evolving some of the relativistic 
components of these complicated scenarios. 

In this large and growing body of work, a significant effort is directed at special 
relativistic magnetohydrodynamics (RMHD). However, with increasing interest in the 
interaction of gravitational and electromagnetic fields, one must also couple the equations 
of relativistic MHD to general relativity either through a curved space background or the 
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dynamical field equations themselves. The earliest attack on this general problem was 
Wilson's pioneering work evolving a rotating, axisymmetric star with a poloidal magnetic 
field | lij. Subsequently, little numerical work was done until very recently with several groups 
developing codes for evolving the general relativistic MHD (GRMHD) equations on fixed 
backgrounds 13]|4||5]|6||7] 111 and in dynamical spacetimes l9lfT0l. 

One of the difficulties alluded to above in the numerical simulation of these sorts 
of systems is that one must solve the GRMHD equations over a large range of time 
and length scales. The computational requirements necessary to adequately resolve 
multiscale phenomena using only a single resolution mesh are often too high for available 
resources. This is particularly true when the full, dynamical GRMHD equations are being 
considered. Adaptive mesh refinement (AMR) therefore becomes crucial in order to reduce 
the computational requirements necessary for modeling such systems. 

There are currently MHD codes in use which incorporate AMR, but most are focused on 
the nonrelativistic problem II llll2lfT3l . Notable among these are the work of Balsara which 
extended a version of constrained transport to AMR. For relativistic MHD, only the work by 
Anninos et al 1 8 1 incorporates adaptive mesh refinement. Theirs is a finite volume approach 
with divergence cleaning. 

This paper describes our algorithm for solving the GRMHD equations with AMR. 
Building on the work presented in 1 141, where the RMHD equations were solved on 
overlapping grids, some key elements of our algorithm are: (1) The Convex Essentially Non- 
Oscillatory (CENO) method for the MHD equations, (2) Hyperbolic divergence cleaning 
for controlling the solenoidal constraint, (3) Berger-Oliger AMR, (4) Weighted Essentially 
Non-Oscillatory (WENO) interpolation for communications from coarse to fine grids, and (5) 
discretization in time via method of lines. 

The CENO scheme is robust, and has three advantages for our code. First, the method 
does not require the spectral decomposition of the Jacobian matrix, making it relatively 
efficient. Central and central-upwind schemes are known to give results nearly identical to 
those of more complicated methods for many problems 11511161 . With the added capability 
of AMR, we find that we are able to efficiently resolve very fine solution features. A second 
advantage of CENO for our purposes is that it is a finite difference, or vertex centered, method. 
As we will solve the Einstein equations with finite differences, using a finite difference fluid 
scheme simplifies coupling the two sets of equations with AMR. The simplification arises 
because fluid and geometric variables are always defined at the same point as grids are refined. 
Finite volume and finite difference grids become staggered with respect to each other as they 
are refined. While we find this simplification advantageous in our present work, we note that 
our AMR code HAD can combine finite difference and finite volume schemes 1171 . A third 
advantage is that ENO schemes are easily extended to higher order accuracy. Our CENO 
scheme can reconstruct fluid variables to both first and second order, resulting in second and 
third order evolution schemes, respectively. 

We choose hyperbolic divergence cleaning 1181 to limit growth in the solenoidal 
constraint on the magnetic field, V • B = 0. Hyperbolic divergence cleaning is simple to 
implement. A single hyperbolic field is added to the system and coupled to the evolution 
equations for B. Divergence cleaning gave good results in earlier tests 1 14 1 and allows us to 
freely choose prolongation methods for AMR. While V • B = is not satisfied to machine 
precision for any particular discrete divergence operator, we find that 1 1 V • B| | does converge 
to zero, and, in the tests presented here, 1 1 V ■ B| | is roughly the same order of magnitude as 
the expected truncation error. 

Mesh refinement is necessary for obtaining accurate numerical solutions for three 
dimensional systems in general relativity, and AMR is essential for complicated problems 
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where the refinement regions can not be guessed a priori. We use the HAD infrastructure 
to provide Berger-Oliger style AMR 1191 . HAD has a modular design allowing one to 
easily implement many different sets of evolution equations, and different modules can be 
combined, for example, to solve both the Einstein and MHD equation simultaneously. HAD 
supports higher order differencing schemes, and our implementation of the MHD equations 
is fully third order accurate 112 01 . Refinement regions can be specified in a variety of ways. 
HAD provides a shadow hierarchy for specifying refinement criteria using truncation error 
estimates, or the user may specify problem specific criteria, such as refining on gradients or 
other solution features. HAD supports different interpolation schemes (we choose WENO 
interpolation for this work) and supports both finite volume and finite difference equations or 
combinations of both. Finally, as discussed below, HAD scales well in strong scaling tests in 
both unigrid and AMR tests. 

Finally, we discretize the continuum equations first in space (creating a semi-discrete 
system) and then discretize in time using the method of lines. This gives us considerable 
flexibility in choosing discretization schemes appropriate for very different types of equations. 
For example, the MHD equations are solved here with high-resolution shock-capturing 
methods, while we might solve the Einstein equations using methods that preserve a 
discrete energy norm. Using the method of lines, we can easily and consistently combine 
these two sets of semi-discrete equations in a uniform time integration. Time integrators 
can be independently chosen for their desired properties or order of integration. For 
example, we choose for this work a third-order Runge-Kutta scheme that preserves the TVD 
condition 1211 . 

The remaining sections of this paper give further details to the algorithm sketched above 
and present code tests in both flat space and on a Schwarzschild black hole background. We 
first present the GRMHD equations used in this work. 



2. The MHD equations in general relativity 

A number of derivations of the GRMHD equations have appeared in the literature, e.g., 1221 
l23M24ll25l l4ll3l fT4l . and thus we simply present the equations to be solved here. The numerical 
methods that we use are very similar to those we have used in our previous MHD work using 
overlapping grids. Thus our presentation here is short, and the reader may refer to that work 
for more detail 1141 . 

The spacetime metric is written in terms of the conventional ADM 3+1 variables, namely 

ds 2 = -a 2 dt 2 + h lJ (dx i + I3 l di)(dx J ' + f3 j dt), (1) 

where a is the lapse, j3 % is the shift, and hij is the 3-metric on the spacelike hypersurfaces. 
Units are chosen such that c = 1 and G = 1. We denote the extrinsic curvature as K a t and 
the Christoffel coefficients with respect to the 3-metric as 3 r^ fc . As our focus in this paper is 
fixed background geometries, we will omit here a discussion of our approach to evolving the 
Einstein equations and only present the relevant matter equations. Future papers will address 
dynamical spacetimes. 

The equations for MHD on a curved background, as in flat spacetime, can be written, for 
the most part, in balance law form, namely 

d t u + d k f k (u) = s(u). (2) 

where it is a state vector, f k are flux functions, and s are source terms. For the current case, 
these are 
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where the quantity (_L T) b is the spatial projection of the stress tensor given in terms of the 
matter fields by 



(±T)\ = v i S b + P.h i b -± 
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In the above, we work from a set of "primitive" variables w = [pQ,v % , P, B^) T consisting 
of the energy density, po, the components of the coordinate velocity of the fluid, v\ the fluid 
pressure, P, and the magnetic field in the fluid frame, BK From these, we define a set of 
conservative variables, u = (D, S b , r, B 3 ') T where the relativistic density D, momentum S b , 
and energy E = r + D are given in terms of the primitive variables by 

D = Wpo, (9) 
S t = [h e W 2 + B 2 ] Vi - (B J Vj) B h (10) 
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and B' J remains unchanged. We have also defined B 2 = B^B 1 , the fluid enthalpy h e = 
po(l + e) + P with e the fluid's internal energy, and the Lorentz factor W = (1 — UjU 1 ) -1 ' 2 . 
Note that spatial indices are lowered and raised by the 3-metric hij and its inverse. Finally, to 
close the system, we assume a T-law equation of state 

P=(T-l)p e, (12) 

where T is the usual adiabatic index. Note that while Eqns. (|3jl-(|6|l are indeed in balance law 
form, Eq. is not. This last equation, of course, is the solenoidal constraint and must be 
dealt with separately. 

One of the benefits of our chosen numerical scheme is that it does not require the full 
spectral decomposition of the system of evolution equations. However, we do find it useful to 
have some information about the possible speeds of some of the waves in our system. This 
information comes by solving for the eigenvalues of the Jacobian matrix, J k ', associated with 
the flux, f k (u), in the k direction where 

df k (u) 
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On doing this for our system for the k direction, say, one gets the general relativistic 
generalization (7) of the seven wave speeds of flat space MHD |26 27 28 29 ]. These include 
the entropy wave and two Alfven waves, 
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and the four ("fast" and "slow") magnetosonic waves, ■ A ± , that are the zeros of the fourth 
order polynomial 
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where c s is the local sound speed and is given by 
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We solve dl6> numerically using the DRTEQ4 routine from the publicly available CERN 
Program Library. The roots from DRTEQ4 are then refined using a Newton-Raphson solver. 

Finally, we note that the MHD equations are written in terms of both the conservative 
and primitive variables. The transformation from conservative variables to primitive variables 
is transcendental, requiring the solution of a single transcendental, nonlinear equation, and is 
outlined in 1141 . 



3. Numerical approach 

This section briefly summarizes the numerical scheme we use to solve the relativistic MHD 
equations. The fluid equations are solved with the Convex Essentially Non-Oscillatory 
(CENO) scheme. CENO is based on a finite difference discretization, which simplifies 
the coupling to the Einstein equations with AMR. Hyperbolic divergence cleaning controls 
growth of error in the solenoidal constraint, and gives some flexibility in choosing other 
components of the numerical algorithm with AMR. 



3.1. CENO 

The CENO scheme was developed by Liu and Osher 11301 to efficiently solve equations in 
balance law form 

d t u + d k f k (u) = «(«), (18) 

where it is a state vector, f k are flux functions, and s source terms. We use the modification 
of CENO for relativistic fluids of Del Zanna and Bucciantini 1311 . Equation dl 81 is solved 
using the method of lines. The semi-discrete form in one dimension is 

^ = - /i+1/ \~ jU/2 +sK), (19) 
at Ax 

where / is a consistent numerical flux. We use both the Lax-Friedrichs flux and the HLL flux 
for the numerical flux. The Lax-Friedrichs flux is 

f F («f+i/ 2 ) + /«i/ 2 ) - «-i/ 2 - u&i/a) 



-;LF 

J i+1/2 



(20) 
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where u^ +1 , 2 and u^_ 1 , 2 are the left and right reconstructed states at x i+ i/ 2 . The HLL 
flux |33 is a central-upwind flux that uses the maximum characteristic velocities for both 
left- and right-moving waves, Xg and A r , respectively, 

A+jV) - Xjf(u r ) + X+Xj(u r - u e ) 



-HLL 



A r A^ 



(21) 



where 

Xj = min(0, X e ) (22) 
A+ = max(0,A r ). (23) 

For highly relativistic flows the Lax-Friedrichs flux gives results very similar to the HLL flux; 
the maximum characteristic velocities approach the speed of light. 

The point-valued fluxes / i+1 / 2 are then converted into consistent numerical fluxes, 

fi+i/2- For a second order scheme, f i+1/2 = /;+i/2> while the correction for the third-order 
scheme is 1211 

/i+1/2 = (1 - -^ (2) ) (24) 

where T>^> is a non-oscillatory second-order difference operator. The operator used in this 
work is specified in 1141 . 

The accuracy of the overall numerical scheme is determined by the reconstruction of the 
fluid states u^ +1 j 2 and uf +1 ^ 2 from the solution known at grid points, i.e., the solution at 
Xi-p, . . . ,Xi, . . . , Xi+ q for integer p and q. Linear and quadratic reconstructions discussed 
below lead to second and third order methods, respectively, for smooth solutions. As is 
commonly done in relativistic fluid dynamics, we reconstruct the primitive variables rather 
than the conservative variables. This is because the conservative to primitive variable 
transformation is transcendental, and computationally rather expensive. 

The reconstruction is performed hierarchically, meaning that a reconstruction of order n 
is calculated from the reconstruction of order n— 1 using a general algorithm. This allows one 
to construct schemes of formally very high order. Thus, we first obtain a linear reconstruction, 
Li, which is then used to create the second order reconstruction. Li is defined on the domain 

[2^-1/2? x i+l/2\ as 

Li(x) = u.i + u[(x - Xi), (25) 
where is the limited slope 

u'i = minmod(Z)_ ui, D + Ui). (26) 

Here we have defined one-sided and centered difference operators as 

n ,Ui±i-Ui u i+1 -u l _ 1 

D±u t = ± , D u, = — , (27) 

Ax 2A.t 

and the minmod limiter is 

!min{aj.} if all > 0, 
max{o)-} ifalla fe <0, (28) 
otherwise. 

The first-order reconstruction, Li(x), is thus equivalent to the linear TVD reconstruction. 

A second order reconstruction proceeds by constructing three candidate quadratic 
functions, Q^(x), which are then compared to Li{x). When the solution is smooth, one of the 
quadratic functions is chosen for the reconstruction. Near discontinuities, however, the linear 
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reconstruction is retained, thus giving solutions similar to TVD schemes for discontinuous 
solutions. The three candidate quadratic functions are 

Qi(x) = Ui +k + D u i+ k(x - x i+ k) + ^D + D-u i+k (x - x l+k ) 2 , (29) 
with k = — 1, 0, 1. Weighted differences with respect to Li(x) are then calculated 

d k (x)=a k (Qfc)-Li(x)). (30) 

The weights a are chosen to bias the reconstruction towards the centered polynomial: 
a = 0.7, and aT 1 = a 1 = 1. When the differences d k all have the same sign, we choose 
the Q k (x) for which d k has the smallest magnitude. When the d k (x) have differing signs, we 
revert to the first order reconstruction. 

Finally, the semi-discrete equations are integrated with the optimal third-order Runge- 
Kutta that preserves the TVD condition 1211 

= u n + AtL(u n ), 

u (3) =lu n + -u^ + -ML{u^), (31) 

4 4 4 v ' 

u n+1 = \u n + + \teL(v?>). 

o o o 

3.2. Hyperbolic Divergence Cleaning 

The time evolution of the magnetic field is governed by (jfji above. However, B is also subject 
to the solenoidal constraint V • B = 0. The continuum evolution equations preserve this 
constraint, although it may be violated in numerical evolutions. These violations can lead 
to unphysical numerical solutions 1331 1341 . Some differencing schemes for the Maxwell 
equations and MHD are designed such that a particular discretization of the solenoidal 
constraint is satisfied to machine precision. These constrained transport methods, naturally, 
do not give solutions that exactly satisfy the continuum constraint, and the magnitude of 
the constraint error can be estimated by using an independent discrete divergence operator. 
Constrained transport methods for classical MHD have been used with AMR [35 36l l37ll38l . 

Divergence cleaning methods are an alternative approach to constrained transport, and 
allow some flexibility in designing the numerical algorithm. Elliptic divergence cleaning 
methods require the solution of a Poisson equation (either explicitly, or implicitly in Fourier 
space), and some common implementations have been reviewed for classical MHD by 
Toth 1 391 an d Balsara and Kim 1401 . Toth reports favorably on divergence cleaning, 
while Balsara and Kim argue that constrained transport performs better for a wider class of 
problems. Hyperbolic divergence cleaning is quite efficient, easy to implement, and usually 
gives good results 1181 . A new field ijj is added to the equations and coupled to the evolution 
equations for B. The field ip acts as a generalized Lagrange multiplier, similar to the A-system 
used in solving the Einstein equations f4Tl . Having some freedom in choosing the equation 
for ip, we choose 

d t B h + d t (B b v l - B l v b ) + h^djtp = 0, (32) 

\d t ip + \ip + V • B =0. (33) 

It can be shown that ip satisfies the telegraph equation, whose solutions are damped, traveling 
waves. The parameters Cf, and c. p control the speed and damping of the constraint propagation, 
respectively. We use Ch — 1 and c p £ [1, 12] in the tests examined here. Using larger values 
of Ch requires smaller Courant factors and did not change results significantly in 1141 . In 
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contrast, the optimal damping factor, c p , is related to the size of the initial shock discontinuity. 
Generally, the larger the shock, the larger the value necessary for c p . Finally, work is 
underway to develop constraint preserving boundary conditions consistent with hyperbolic 
divergence cleaning for the MHD equations 1421 . 

3.3. Adaptive Mesh Refinement 

AMR provides the ability to add grid refinement where and when needed. This need is 
determined by some refinement criterion. At any given level of resolution, points which 
meet this criterion are flagged and a new, finer level is created which includes all such flagged 
points. Similarly, when no points are flagged, the level is removed. For the tests presented 
here, the maximum number of levels of refinement was limited to two. In other words, our 
runs have grids with three different resolutions. 

The fine and coarse grids communicate in AMR through prolongation and restriction. 
Fine grids are created by interpolating the solution from a parent grid (prolongation), 
and the fine grid solution is communicated to coarser grids through restriction. In 
prolongation we interpolate the conservative variables onto finer grids using third-order 
WENO interpolation 1431 1141 . This interpolation scheme is designed to work well with 
discontinuous functions by adjusting the interpolation stencil to the local smoothness of the 
function. This avoids oscillations near discontinuities, which often cause primitive variable 
solvers for relativistic fluids to fail. For restriction on vertex centered grids, the fine grid 
values are copied directly to the coarse grid (direct injection). If a point on a vertex centered 
coarse grid is also found on two or more finer grids, the restriction operation averages the 
values on the finer grids for the solution at the coarser grid point. 

The distributed AMR infrastructure that we employ is HAD. HAD is a F77 based Berger- 
Oliger 1 1 ') i type AMR package presented in |44j using the message passing interface (MPI) 
for distributed parallelism. HAD has a modular design, allowing one to solve different 
sets of equations with the same computational infrastructure. Unlike many other publicly 
available AMR toolkits, including E3 |46j |47] EHJ |49| , HAD is natively vertex centered. This 
can be advantageous in numerical relativity, as many difference schemes for the Einstein 
equations are vertex centered. Support for cell centered grids in HAD is also available. 
HAD supports subcycling of grids in time for full space-time AMR and can in principle 
accommodate arbitrary orders of accuracy in both space and time. An example has been 
shown using third order accurate AMR simulations 1201 . The HAD clustering algorithm is 
Berger-Rigoutsos 1501 : the load balancing algorithm is the least loaded scheme 1511 . 

Considerable flexibility is provided in developing a refinement criterion for the HAD 
infrastructure, and many were explored in conjunction with the numerical tests presented 
here, such as refining on gradients in the density, pressure or the magnetic field. All criteria 
resulted in similar adaptive mesh hierarchies for the relativistic rotor and spherical blast wave 
tests; consequently, those tests center the refinement on the evolving shock front. A shadow 
hierarchy 1 52 1 has recently been added to HAD for truncation error estimation and was used as 
the refinement criterion in the spherical accretion tests of a fluid falling onto a Schwarzschild 
black hole. 

4. Numerical Results 

In this section we examine three relativistic MHD test problems and one accretion test 
problem using AMR. The MHD test problems are selected because their solutions have very 
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Balsara blast wave parameters 
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Table 1. The initial parameters for the Balsara blast wave. The discontinuity is initially placed 
at x = and the system is evolved along the x axis. The adiabatic index is V = 5/3. 
The exact solution to this problem is given in 1531 . The Courant factor used in the numerical 
solution is 0.2. 



sharp features that require high resolutions. The accretion test problem demonstrates the 
AMR capabilities with a curved space background and an excision region. 

The first problem is a one-dimensional blast wave problem introduced by Balsara II II . 
We compare the unigrid shocktube results with the exact solution, and then compare unigrid 
results with the resolution-equivalent AMR results. The second and third problems are three- 
dimensional extensions of standard two-dimensional tests: a spherical blast wave and a 
spherical relativistic rotor 1281 1101 fl4l . We present unigrid and AMR results for these test 
problems, and discuss the effects of hyperbolic divergence cleaning. The last problem is the 
accretion of a fluid onto a Schwarzschild black hole. We numerically recover the steady state 
solution using AMR with a shadow hierarchy as the refinement criterion. Finally, we present 
scaling results for parallel unigrid and AMR runs. In the tests below we use the discrete Li 
norm 
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where u is a discrete function defined at N locations, 



(34) 



4.1. Riemann problem test 

Balsara introduced several test Riemann problems for relativistic MHD II II . and we choose 
here his third blast wave problem for its very narrow features as a test of our AMR. The initial 
parameters for this Riemann problem are given in tabled The exact solution for this problem 
is given by Giacomazzo and Rezzolla 1531 . and is plotted in the figures below for comparison. 

The blast wave problem is implemented using the three-dimensional HAD infrastructure 
for AMR and simulated along a coordinate axis. Figure^shows the blast wave with a strong 
initial pressure difference centered at x = evolved along the x-axis at several unigrid 
resolutions. The plots show po, v x , v v , and B v at time t = 0.4. The unigrid simulations 
were performed in the x direction on a domain of [—0.5, 0.5] with a Courant factor of 0.2. 

A series of two-level AMR simulations of the blast wave test were conducted with 
refinement centered on the shock propagation. Figure |2] compares a unigrid simulation 
with a two-level AMR simulation where the resolution of the finest mesh in the AMR 
hierarchy is the same as that in the unigrid simulation. The AMR capability to reproduce 
the resolution-equivalent unigrid result depends on how well the refinement region tracks the 
shock throughout the time history of the evolution. In our tests we observe the AMR blast 
wave simulations reproducing the resolution-equivalent unigrid result to within 0.1%. 
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Figure 1. Unigrid simulation results for the relativistic MHD Balsara blast wave test at time 
0.4 showing po,v x ,v y , and B y . The z components of B and v are identical to their respective 
y components. The simulations were performed along the x axis using four resolutions on a 
domain of [—0.5,0.5]. The base resolution was ho = 6.25 X 10 — 4 . The exact solution to 
this problem is found in 1531 . This problem is an excellent candidate for AMR because of the 
high resolutions required to adequately resolve the different waves. 
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Figure 2. Adaptive mesh refinement results for the Balsara blast wave. Figure|2ja) shows the 
difference in po between the finest mesh in a two-level AMR simulation and the equivalent 
unigrid simulation. The finest mesh of the AMR simulation has the same resolution as the 
unigrid case, h = 1.25 X 10 -4 . Simulations were performed on a domain of [—0.8, 1.2]. 
The unigrid required 16000 points while the AMR required 5423 points. For reference, po antl 
the layout of the AMR mesh hierarchy at t = 0.4 are shown in the window inset. Figurel2*tbi 
compares two solutions of roughly equal computational cost: one unigrid using 6250 points 
with resolution h = 3.2 X 10~ 4 and the other AMR using 5423 points with maximum 
resolution h = 1.25 X 10~ 4 . The AMR simulation produces significantly more accurate 
results for po than the unigrid solution at the same computational cost. 



4.2. Spherical blast wave 

The spherical blast wave consists of a uniform fluid background with a small spherical region 
where the pressure is 10 6 times larger than the background. The background pressure is 10~ 2 , 
and P = 10 4 inside a central sphere of radius 0.08. This is the three-dimensional extension of 
the cylindrical blast wave studied in 1281 1 101 PHI . The parameters for the spherical blast wave 
are given in table|2] 

We first calculate the solution on a single uniform grid, and then draw comparisons 
to the AMR results. Figure [3] shows z — cuts of the uniform grid solution at t = 0.4, 
and figure |4] shows line plots of the pressure, P, along the x- and y-axes at three different 
resolutions. Adaptive mesh refinement substantially improves performance for the spherical 
blast wave while returning results nearly identical to the unigrid result. Figure [5] shows the 
resulting mesh hierarchy and pressure at time t = 0.4 for the spherical blast wave in an AMR 
simulation with two levels of refinement. This AMR simulation uses hyperbolic divergence 
cleaning and is the AMR equivalent of the divergence cleaning unigrid simulation in figure|3] 
The refinement criteria are set to center refinement on the shock. The relative simplicity 
of this solution — the outgoing shock is the dominant feature of the solution-allows many 
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Table 2. The initial parameters for the spherical blast wave. The data consists of a uniform 
fluid background with the pressure set to 10 4 inside a sphere of radius 0.08 centered at the 
origin. The adiabatic index is T = 4/3. The domain of simulation is {x, y, z} 6 [—1, 1]. 
The Courant factor is 0.3. 



With Cleaning Difference 




P 9.07e-03 3.27e+01 AP -1.92e-01 3.48e-01 



Figure 3. Unigrid simulation results for the spherical blast wave at t = 0.4, showing a slice 
along the 2 = plane. The x axis is the horizontal direction. The pressure found using 
hyperbolic divergence cleaning is shown on the left. The difference between the pressure 
found with and without hyperbolic divergence cleaning is shown on the right. This gives an 
estimation of the relative errors that arise in free evolutions. The simulations were performed 
using a resolution of h = 0.006410 on a domain of {x, y, z} £ [—1, 1]. 

refinement criteria to produce similar mesh hierarchies. The AMR simulation in figure[5]was 
performed on 32 processors and was a factor of eight times faster than the equivalent unigrid 
simulation. Like the Balsara blast wave case examined in l4.ll adaptivity significantly reduces 
the computational overhead required to adequately resolve the multiscale features that appear 
in the simulation. 

Finally, we monitor the violations of the solenoidal constraint during both free evolutions 
and evolutions with hyperbolic divergence cleaning. Figure[6]shows the magnitude of the L2 
norm of V • B at three different resolutions with and without divergence cleaning. There 
are some subtleties in interpreting L2 norms of the constraint violation, which arise primarily 
because Richardson-like convergence can not be defined for discontinuous functions 1141 . 
However, it appears here that the constraint violations are propagated at roughly the same 
velocity as the out-going shock, and the difference between free and cleaned evolutions is 
small. 
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(a) 



(b) 



Figure 4. The pressure for the spherical blast wave at t = 0.4 shown at three resolutions along 
the a;-axis (left frame) and the y-axis (right frame). The window insets provide a closer view 
of the shock front amplitudes at the different resolutions. The base resolution is ho = 0.0128. 



4.3. Relativistic rotor 

The relativistic rotor test case starts with a rigidly rotating fluid and evolves it in the presence 
of a magnetic field. This problem is discussed and examined in 2 + 1 dimensions in 11141 1281 - 
Here we examine the relativistic rotor in 3 + 1 dimensions, confining the initially rigidly 
rotating fluid to a sphere of radius 0.1 with the angular momentum vector pointing in the +z 
direction. The fluid is initially rotating with an angular velocity of 9.95. The initial data and 
relevant evolution parameters are given in table [3] 

Results using a uniform computational grid form the standard against which we measure 
the AMR results. The L2 norms of the V ■ B constraint violation as a function of time using 
three different unigrid resolutions are shown in figure0 comparing results obtained both with 
and without hyperbolic divergence cleaning. Hyperbolic divergence cleaning significantly 
improves the constraint preservation in the relativistic rotor case. 2-D slices of the pressure 
along the z = plane at time t = 0.4 are shown in figure[S] 

Adaptive mesh refinement results are presented in figures |9] and [TO] These figures 
present a two-level AMR simulation with refinement centered on the shock front. This 
AMR simulation required five times fewer CPU hours than the equivalent unigrid simulation. 
Figure|9]shows the resulting pressure and mesh hierarchy at time t = 0.4. FigureflQlcompares 
the difference along the x and y axes of the pressure between the AMR and equivalent unigrid 
simulation. The AMR simulation was found to reproduce the unigrid results to well within 
0.1%. 
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Figure 5. A wireframe slice along the z = plane of P for the spherical blast wave using 
hyperbolic divergence cleaning and two levels of adaptive mesh refinement. In each plot, the 
x axis is the horizontal direction. The top frame illustrates the domain decomposition of the 
system across 32 processors. Only two separate resolutions are distinguishable, but a third 
resolution becomes apparent in close-up, shown in the the bottom frame. The simulation 
reproduces the peak unigrid amplitude to within 0.4%. This AMR simulation is eight times 
faster than the equivalent unigrid simulation (shown in figure |3j and was performed with 
a maximum resolution of h = 0.006410 on a domain of {x,y,z} £ [—1,1] using 32 
processors. 
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Figure 6. The L2 norm of V ■ B as a function of time for the spherical blast wave at three 
resolutions for both free evolutions and evolutions with hyperbolic divergence cleaning. The 
base resolution is ho = 0.025641. The divergence cleaning parameters are c/j = 1 and 
c p = 12. 
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Table 3. The initial parameters for the relativistic rotor. The data consists of an initially rigidly 
rotating fluid inside a sphere of radius 0.1 centered at the origin with a magnetic field. The 
fluid is rotating with ui = 9.95 around the z axis. The adiabatic index is T = 5/3. The 
domain of simulation is [—1, 1] in each of the x,y, and z directions. The Courant factor used 
in the presented relativistic rotor simulations was 0.2. 
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Figure 7. The L2 norm of the V ■ B = constraint violation as a function of time for 
the relativistic rotor at three resolutions comparing results obtained with free evolution and 
hyperbolic divergence cleaning. Hyperbolic divergence cleaning has a profound impact on 
constraint control. The base resolution is ho = 0.025. The hyperbolic divergence cleaning 
parameters used were = 1 and c p = 1, 



4.4. Accretion onto a Black Hole 

Numerical simulation of fluid accretion onto a Schwarzschild black hole provides a test of the 
AMR infrastructure using a curved space background and an excision region. The steady state 
solution is given by Michel |54| and has been explored previously using ingoing Eddington- 
Finkelstein coordinates 1551 . This solution describes the continuous spherical accretion of a 
fluid onto a black hole. 

We set a fixed black hole metric using ingoing Eddington-Finkelstein coordinates. To 
avoid the singularity inside the black hole, we implement a cubic excision region. The 
excision cube is located at the center of the grid and has a half width of 0.3M. The boundary 
condition at the excision region is a copy condition; points next to the excision region are 
simply copied into the excision region when necessary for reconstruction. The mass of the 
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4.14e-01 1.67e+00 Ap -7.37e-02 4.86e-02 

Figure 8. Unigrid simulation results for the relativistic rotor at time 0.4, showing a slice along 
the 2 = plane. The x axis is the horizontal direction. The pressure found using hyperbolic 
divergence cleaning is shown on the left. The difference between the pressure found with and 
without hyperbolic divergence cleaning is shown on the right. This gives an estimation of the 
relative errors that arise in free evolutions. The simulations were performed using a resolution 
of h = 0.00625 on a domain of {x, y, z} G [-1, 1]. 




Figure 9. A wireframe slice along the z = plane of the pressure for the relativistic rotor 
using hyperbolic divergence cleaning and two levels of adaptive mesh refinement. The x 
axis is the horizontal direction. This simulation reproduces the unigrid solution to better 
than 0.1%. See figure fTol This AMR case required five times fewer CPU hours than the 
comparable unigrid case (shown in figure |5J and was performed with a maximum resolution 
of h = 0.00625 on a domain of {x, y, z} € [—1, 1] using 32 processors. 
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Figure 10. The absolute difference in pressure between a two level AMR simulation and the 
equivalent unigrid simulation for the relativistic rotor case at time 0.4. Figure [Tol a^ plots the 
difference in pressure between the AMR and unigrid results along the x axis; figure Hol b) plots 
the difference along the y axis. The AMR simulation is the same as that shown in figure l9l 
The unigrid simulation is the same as that shown in figure[8] The finest resolution meshes of 
the two level AMR system had a resolution of h = 0.00625, equivalent to the unigrid mesh 
resolution. The results are identical to better than 0.1%. For reference, the pressure along the 
x and y axes are plotted in the window insets. 



black hole, M, is set to one and the black hole is placed at the center of the grid. The sonic 
radius, r c , is selected to be 400M with a density p c — 0.01. The domain of simulation is 
{x, y, z} G [-15M, 15M]. The Michel steady state solution is found following the procedure 
described in 1551 and the outer boundary is kept fixed at this solution, providing a continual 
influx of mass onto the black hole. For radius r > 2.5M the Michel steady state solution is 
set as initial data; for r < 2.571/ the initial data are set to be po = 0.1, P = 0.1, and v l = 0. 
The fluid falls onto the black hole and eventually reaches steady state. A comparison of the 
Michel steady state solution and numerical solution at t = 50M is given in figure ^2 The 
AMR grid structures at times t = 0M and t = 50 M are given in figure [21 A convergence 
test for po is presented in figure [T3l 

4.5. Scaling 

Unigrid and mesh refinement parallel scaling tests for the spherical blast wave are given in 
figureO The results presented are the strong scaling results; the global problem size was kept 
constant while the number of processors varied. Strong scaling tests are problem dependent 
and vary according to the size of the global problem selected for investigation. However, they 
also give the most direct indication of performance speed-up across a wide range of processors 
for a particular problem. Weak scaling tests are where the global problem size is increased as 



Relativistic MHD with Adaptive Mesh Refinement 



19 




Figure 11. This figure compares po from the numerical steady state solution of accretion onto 
a black hole with the Michel solution along the x— axis of the computational domain. The fluid 
is initially set to the Michel solution for radius r > 2.5M and constant pressure and density 
for r < 2.5 A/. The system is then evolved until steady state is reached. The numerical result 
shown here is at t = 50M with a finest resolution of 0.075M . The excision region of the 
black hole is in the center of the grid. The AMR grid structure of the simulations is shown in 
figure H2l 




t = OM t = 50M 

Figure 12. This figure shows po and the AMR grid structures at t = OM and t = 50M 
along the x — y plane. The refinement criteria is the shadow hierarchy for truncation error 
estimation. The fluid is initially set to the Michel solution for radius r > 2.5M and constant 
pressure and density for r < 2.5M. The system is then evolved until steady state is reached. 
The cubical excision region is highlighted in the center of the grid on the left. 
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Figure 13. This figure plots log(Apo/po) along the x-axis of the computational domain 
where Apo is the difference in po between the numerical steady state accretion solution and 
Michel solution at t = 50M. The convergence test consists of successively enabling another 
AMR refinement level to get a finer resolution. The resolution reported in the legend of the plot 
is the finest resolution present in the simulation. Only for the lowest resolution case, ht\, is 
the resolution constant across the entire grid. Cases hf2 and hf^ contain multiple resolutions 
across the domain. At locations where simulations share the same resolution, they also display 
the same error, modulo AMR boundary effects. From x € [2.1 M, 15 M], simulations hf2 
and hfs have the same resolution. From x S [3.9 M, 15 M] all three simulations share the 
same resolution. The coarsest resolution for these simulations is ho = 0.3 M. 



the number of processors increases in order to keep the problem size local to each processor 
constant. Weak scaling tests were also performed on the spherical blast wave problem on 
16-256 processors without showing any significant performance degradation. 

In the unigrid strong scaling tests of figure \\4\ a) a 121 3 spherical blast wave problem 
was evolved for 80 iterations on 1-128 processors. Speed-up is defined as 

Run time on one processor 

spccdup(n) = — ; . 

Run time on n processors 

As the number of processors increases, the communication eventually overshadows the local 
processor computation. For the test problem size examined, this begins on > 64 processors. 
For comparison, strong scaling results for a different unigrid TVD MHD code are given in 
1561 . where communication overhead saturation occurs on > 64 processors using a 240 3 
mesh size and 8 processors the base scaling value. 

In the mesh refinement strong scaling tests of figure fl4T b) a 81 3 spherical blast wave 
with a single level of mesh refinement was evolved for 30 iterations. The base number of 
processors for scaling measurement was 8 on account of memory considerations. Results are 
shown for 8-80 processors. 

5. Conclusion 

We have presented three flat space relativistic MHD tests and one fluid accretion test using 
vertex-centered distributed adaptive mesh refinement and the approximate Riemann solver 
algorithm for the GRMHD equations presented in 1141 . Each of the three relativistic MHD 
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Figure 14. This figure shows strong scaling results for single grid and mesh refinement. The 
spherical blast wave initial data were run for a fixed problem size as the number of processors 
is varied. The left frame shows the unigrid strong scaling results, and the right frame includes 
mesh refinement. For the unigrid scaling, the data were evolved for 80 iterations, and the global 
grid size was 121 3 . For this problem size, the communication overhead begins to overshadow 
the local process computation on > 64 processors. For the mesh refinement scaling, thirty 
iterations were performed on a coarse grid of size 81 3 and a single level of refinement. Since 
the test problem would not fit in memory on a single processor, speed-up was measured using 8 
processors as the base value. All tests were performed on an Intel Pentium IV 3.0 GHz cluster 
with Myrinet. 

tests, including the Balsara black wave and the spherical shock and relativistic rotor, have 
sharp features requiring high resolutions. In each of these cases, substantial performance 
gains of AMR versus unigrid were observed. Two level AMR simulations required between 
5-8 times fewer CPU hours than the equivalent unigrid cases. AMR results reproduced the 
unigrid results, often better than 0.1%. 

Hyperbolic divergence cleaning was examined in connection with the spherical blast 
wave and relativistic rotor cases. It had a positive impact on constraint control in both cases. 
The impact was especially pronounced in the relativistic rotor case, which has more features 
interior to the outgoing wave front than the spherical blast wave. Hyperbolic divergence 
cleaning worked equally well in both unigrid and AMR tests. 

Fluid accretion onto a Schwarzschild black hole tests the code using a curved space 
background and a black hole excision region as well as using AMR. In this test, a shadow 
hierarchy for truncation error estimation was used as the AMR refinement criterion in 
recovering the known steady state solution. 

Parallel performance measures were presented in connection with the spherical blast 
wave. Speed-ups achieved using the spherical blast wave were reported both with unigrid 
and mesh refinement simulations. Performance speed-ups were found on up to at least 128 
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processors with unigrid and up to at least 80 processors with mesh refinement. By its very 
nature, the strong scaling test is problem dependent: the problem size is fixed while the 
number of processors is varied. This is in contrast to the weak scaling tests often presented 
in numerical relativity, where the problem size per processor is fixed. Strong scaling tests, 
however, address the real-world questions of how long it takes to solve a particular problem, 
and how to do it most efficiently. 

Having presented these tests, we now turn to some questions of astrophysical interest 
mentioned in the introduction. In particular, we have included fully dynamical general 
relativity in our code using the Einstein equations specified in 1571 . In future work we hope 
to present evolutions of TOV stars as well as rotating, magnetized neutron stars. Recent work 
suggests interesting effects of general relativity with rotation on supermassive polytropes and 
the bar mode instability |58|. The addition of magnetic fields to these systems may suggest 
new questions and provide new insight into magnetic astrophysical phenomena. One part of 
this question includes understanding not only the interior of a magnetized rotating neutron 
star, but its magnetosphere as well. Ideal MHD codes based on Godunov-type schemes 
frequently encounter difficulties when B 2 ^> po, as relatively small truncation errors in 
the evolution of the conserved variables lead to large fractional errors in computing the 
internal energy density and other primitive variables. While we have tried to create a robust 
primitive variable solver, this difficulty can not be avoided for high-resolution shock-capturing 
schemes. Therefore, a full numerical study of such a star and magnetosphere may require 
coupling the equations of ideal MHD for the interior solution with the equations of force-free 
electrodynamics for the exterior. We are actively pursuing this question. 

Another area to be targeted in future papers is constraint preserving boundary conditions 
for MHD. Currently we use the conventional outflow boundary conditions. This could be 
improved by using outer boundary conditions that are constraint preserving. Additionally, 
for some systems we wish to require that no incoming modes enter the domain. To this end 
it would be useful to construct the full spectral decomposition of our system in order to be 
able to determine ingoing and outgoing modes. Work in this direction has already begun and 
shows promise. 
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